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Abstract. A model of the quiet time middle magnetotail is developed using a consistent 
orbit tracing technique. The momentum equation is used to calculate geocentric solar 
magnetospheric components of the particle and electromagnetic forces throughout the 
current sheet. Ions generate the dominant x and z force components. Electron and ion 
forces almost cancel in the y direction because the two species drift earthward at compara- 
ble speeds. The force viewpoint is applied to a study of some substorm processes. Genera- 
tion of the rapid flows seen during substorm injection and bursty bulk flow events implies 
substantial force imbalances. The formation of a substorm diversion loop is one cause of 
changes in the magnetic field and therefore in the electromagnetic force. It is found that 
larger forces are produced when the cross-tail current is diverted to the ionosphere than 
would be produced if the entire tail current system simply decreased. Plasma is acceler- 
ated while the forces are unbalanced resulting in field lines within a diversion loop becom- 
ing more dipolar. Field lines become more stretched and the plasma sheet becomes thinner 
outside a diversion loop. Mechanisms that require thin current sheets to produce current 
disruption then can create additional diversion loops in the newly thinned regions. This 
process may be important during multiple expansion substorms and in differentiating 
pseudoexpansions from full substorms. It is found that the tail field model used here can 
be generated by a variety of particle distribution functions. However, for a given energy 
distribution the mixture of particle mirror or reflection points is constrained by the consis- 
tency requirement. The study of uniqueness also leads to the development of a technique 
to select guiding center electrons that will produce charge neutrality all along a flux tube 
containing nonguiding center ions without the imposition of a parallel electric field. 

length of 0.7 R £ in this region. All the numerical calculations were 
fully three-dimensional (3-D). However, since our standard mag- 
netic field model has significant y dependence only near and sun- 
ward of x = -10 /?£, the model is essentially 2-D in the region of 
interest. 

Since the COT analysis finds particles that create a model mag- 
netic field, it is important to start with a model that is representative 
of the actual magnetotail. T89 is an empirical model involving 
many parameters that were adjusted to fit a large collection of sat- 
ellite data. The data were sorted according to geomagnetic indices 
such as Kp but not according to substorm phase. As a result, the 
model tail represents an average over both the thin stretched con- 
figuration that is seen before a local substorm expansion and the 
dipolar configuration seen after expansion. It therefore is not nec- 
essary that the model will provide a good representation of the 
Earth's field at any one time. 

One important question is whether the model is consistent with 
steady convection. Erickson and Wb//[1980] pointed out that adia- 
batic fluid convection in common magnetotail models will result in 
plasma pressures in the inner plasma sheet that are much higher 
than is observed. Some models have specifically been developed to 
be consistent with quasi static convection [ Hau , 1991]. This shows 
that it is possible to develop a steady state magnetic field model in 
which forces are nearly balanced. Chen and Wolf[ 1 993] developed 
a nonuniform model containing small plasma bubbles that are 
accelerated primarily in the distant tail and decelerated primarily in 
the near-Earth plasma sheet. This model was constructed so that it 
can explain the observed pressure profile. Finally, it is possible that 
the plasma sheet is usually evolving smoothly as in many MHD 


1. Introduction 

1.1. Consistent Orbit Tracing Method 

The consistent orbit tracing (COT) technique [Larson and Kauf- 
mann , 1996] (hereafter referred to as LK96) was used to construct 
a steady state model current sheet. In a consistent model the ions 
and electrons carry the electric current needed to produce the mag- 
netic field in which their orbits were traced. The COT method 
begins with a preselected magnetic field and a search is made to 
find particles that will produce this field. In contrast, most simula- 
tions begin with preselected particle sources and boundary condi- 
tions and calculate the resulting fields. 

The standard magnetic field plus uniform cross-tail electric field 
model described in LK96 was used here. The magnetic field was 
adjusted so that the adiabaticity parameter defined by k = 
/? min /p max [Buchner and Zelenyi , 1989] and the total cross-tail 
sheet current density K y (x) in A/m, which is the volume current 
density j (x, z) in A/m 2 integrated through the thickness of the 
current sheet, were similar to those at y = 0 in the Kp = 4 version of 
the Tsyganenko [1989] or T89 model. In the above definition /? mjn 
is the minimum radius of curvature of a magnetic field line and 
p max is the maximum gyroradius of an ion. Both /? mjn and p max 
are found at z = 0. The region studied in detail was —20 Rg < x < 
-14 R e , 0 < Izl < 2 R e in geocentric solar magnetospheric coordi- 
nates. The standard model current sheet has a characteristic z scale 
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models [Walker et al ., 1993; Birn et al., 1996]. In this scenario the 
field approaches a quasi-static configuration [ Hau , 1991] only dur- 
ing periods of steady convection [Sergeev et al. % 1996a]. These top- 
ics will be discussed in more detail in a paper that is devoted to a 
study of the energy equation, which describes the evolution of the 
pressure tensor as plasma flows toward the Earth. 

Groups of 1000 ions each were traced in the model fields to 
carry out a COT analysis. A starting energy, pitch angle, and phase 
angle were randomly selected for each ion using a 5-keV Max- 
wellian parent population. The starting x and z positions were ran- 
domly distributed over one of the spatial boxes. Each group used in 
LK96 also was composed of 1000 ions, but those ions were 
selected from various monoenergetic parent distribution functions. 
No explicit boundary conditions were imposed at the edges of the 
region of interest. Particles drifted from a physical source some- 
where tailward of x = -20 R E to a sink somewhere earthward of x = 
-14 /?£. Orbits were traced both forward and backward in time 
until the particle drifted so far it could never return to the region of 
interest. 

The energy of each ion changed as the particle drifted and dif- 
ferent groups were started at different locations. The final energy 
distribution in a given spatial box therefore was broader than the 
starting distribution. Full 3-D velocity distribution functions were 
retained in 1 20 spatial boxes for each group of ions that was traced. 
The boxes were 0. 1 -R E wide in the z direction and l-/? £ wide in the 
x direction. The model is symmetric about z = 0. Only one very 
wide box was used in the y direction because no useful information 
would be gained by keeping multiple y boxes in a region that has 
essentially no y dependence. 

Electron currents were added to the current calculated for each 
group of ions. A number of ion-plus-electron groups then were 
combined to form a plasma that approximately carries the current 
which is needed to produce the magnetic field. The combination of 
groups was done with a least squares fitting method that considered 
only the y component of j. 

The decision to use 1000 ions for each group was made by sim- 
ply increasing the number of particles per group until repeating the 
orbit tracing process with multiple groups which were randomly 
selected from the same parent population produced almost identi- 
cal current patterns. Since only j y was used for the fitting routine, 
there was no advantage in increasing the size of each group at this 
point in the analysis. Larger groups were used for some studies, 
such as those which needed accurate values of off diagonal ele- 
ments in the pressure tensor. 

A total of 20 to 40 groups of ions usually were traced to carry 
out a COT analysis. However, an average of only about eight 
groups was kept to produce the final plasma. Keeping more groups 
produced better agreement between the j y carried by particles and 
the j that is needed to produce the model field. However, our 
principal goals involve learning more about the physical structure 
of the current sheet rather than producing the closest possible fit to 
a preselected model. We therefore required that each group in the 
final plasma must be needed at a 95% confidence level to improve 
the agreement between the model j y and the desired j y . 

1.2. Background, Goals, and Overview 

The advantages of the COT method are it produces approxi- 
mately consistent kinetic current sheet models which satisfy the 
Vlasov equation; the use of very accurate orbit tracing assures that 
all nonguiding center aspects of particle trajectories are included; 
and 3-D ion velocity distribution functions ( r, v) are generated 
for each of the spatial boxes. The availability of f. (r, v) permits a 


detailed study of the structure of the neutral sheet region, where ion 
orbits are very complex. It is easy to evaluate all the fluid parame- 
ters and derivatives that are needed to calculate forces and energy 
flows and to make comparisons with observations. 

The most important disadvantage of the COT technique is that it 
only generates a steady state model. The topics of stability, fluctua- 
tions, and time evolution are not addressed by these calculations. 
The COT results may provide good starting distributions for future 
kinetic simulation and instability studies. For practical reasons the 
present COT results covered a relatively small region of the mag- 
netotail. The results in this paper therefore cannot be used to locate 
the solar wind and ionospheric sources of plasma sheet particles, as 
was done using the large-scale kinetic orbit tracing calculations 
[Peroomian and Ashour-Ahdalla , 1995], The limited region that 
was modeled using the COT technique also means that it is not 
possible to integrate throughout a flux tube from the ionosphere to 
the equator. This limitation makes it hard to compare the present 
results to bounce averaged calculations. 

Methods that are similar to the COT technique were used by 
Eastwood 11975], Other early orbit tracing studies relevant to the 
COT technique include Speiser [1965, 1970], Alekseyev and Kro- 
potkin [1970], Sonnerup [1971], Francfort and Pellat [1976], 
Cowley and Pellat [1979], Wagner et al. [1979], and Gray and Lee 
[1982]. The application of techniques from chaos theory by Buch- 
ner and Zelenyi [ 1 986] and Chen and Palmadesso [1986] were par- 
ticularly important. These and other early studies are described by 
Kaufmann and Lu [1993] and Kaufmann et al. [1994] where a ver- 
sion of the COT method suitable for 1-D magnetotail models was 
developed. The I-D studies found that it is easy to generate COT 
current sheets that produce a modified Harris [1962] model mag- 
netic field. However, the resulting current sheets were not consis- 
tent with observations by satellites in the magnetotail. The 2-D 
models used in LK96 produced current sheets that agreed with 
observations. 

The first goal of the present paper is to see how accurately 
forces are balanced in the COT model and to calculate the relative 
importance of each term in the momentum equation. Average mag- 
nitudes of these terms were estimated long ago using quiet time 
measured plasma parameters and fields [Cole and Schindler , 1972; 
Rich et al. , 1972]. The present paper gives a more detailed compar- 
ison using the full pressure tensor to calculate the x and z spatial 
dependence of each Cartesian component of the force terms. It also 
is important to check the reliability of the COT results that use 
derivatives of the pressure tensor. The accuracy of such derivatives 
must be assessed for a future study of energy flow, heat flow, 
entropy changes, and equations of state. Section 2 is devoted to a 
study of the momentum equation and an evaluation of terms 
derived by the COT analysis. 

Events that occur near substorm onset are studied in section 3 as 
a demonstration of the usefulness of the force viewpoint. Substorm 
processes in the magnetotail usually are described in terms of the 
disruption of current in the near-Earth plasma sheet and the diver- 
sion of current to the ionosphere. The alternative force viewpoint 
emphasizes the unbalanced forces that develop near substorm 
onset. The fact that changes in both the local electromagnetic force 
and the local bulk flow velocity can be measured in data from a 
single satellite also may make it possible to directly observe some 
consequences of force imbalance. It is much more difficult to 
directly measure changes in a full substorm current system. 

The above results are summarized and compared with recent 
orbit tracing and kinetic analyses in section 4. The appendices 
address some questions of uniqueness. It is seen that a range of dis- 
tribution functions could produce the preselected 2-D magnetic 
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field. The more complex problem involving 3-D fields with strong 
y dependence is not treated. It also is shown that the method used 
to add electrons and to calculate the associated parallel electric 
field is not unique. A technique is derived to adjust the equato- 
rial guiding center electron pitch angle distribution to assure charge 
neutrality with a given nonguiding center ion density distribution 
and £„ = 0. 

A companion paper [Kaufmann et al ., this issue] considers non- 
guiding center aspects of the COT model. This companion study 
emphasizes the importance of thin current sheets and weak mag- 
netic fields near z = 0. Nonguiding center mechanisms that may 
disrupt cross-tail current in a thin current sheet without signifi- 
cantly decreasing the tail plasma population also are described. 
Future papers will be devoted to an analysis of the energy equation 
and a comparison of the calculated 3-D velocity distribution func- 
tions with observations. 

2. Force Balance 


responding to the end of the time step. Velocity space was divided 
into 30 boxes* 15 for each positive and 15 for each negative Carte- 
sian component of velocity. The result is a 27,000-point 3-D veloc- 
ity distribution function for each of the 120 spatial boxes. 

The stress tensor T, contains everything that is needed concern- 
ing the kinetic energy of plasma particles. The separation on the 
right side of (1) into , which contains information about the 
thermal kinetic energy, and the Reynolds stress tensor, which con- 
tains bulk flow information, will be used throughout this work. 

In contrast, the electromagnetic stress tensor depends only upon 
the model electric and magnetic fields f Jackson , 1975, p. 239] 


rp EM 


%E a E^^B a B p-t 


e„EE + f i-B-B 




(4) 


Forces produced by the particles whose orbits were traced to 
generate the model current sheet will be compared with forces 
exerted by the preselected fields. 


2.1. Stress Tensors 

Forces can be calculated by taking derivatives of various stress 
tensors. The particle stress or kinetic tensor for particles of species 
s is [ Spitzer , 1962; Rossi and Olbert , 1970] 

T S, a (5 = m sj v a v p fs ( r - v ) d ~ v = p s% + P A <A p (1 > 

where 

V i,a = J.f v a/ S ( r ’ v )^ v < 2 > 


2.2. Force Equations 

A model current sheet generated by the COT technique ideally 
results in force balance. This can be seen by noting that the Vlasov 
equation is based on the assumption that each particle follows an 
unperturbed orbit in the average B and E fields produced by all 
other ions and electrons. The COT technique explicitly follows 
such unperturbed orbits. This method therefore is equivalent to 
solving the Vlasov equation provided the final plasma sheet is con- 
sistent in the sense described above. The momentum equation is 
derived by taking the first velocity moment of the Vlasov equation. 
For nonrelativistic particles with charge q and neglecting gravity 
this gives [ Spitzer , 1962; Rossi and Olbert , 1970] 


is the bulk flow velocity of the species 

(vp-V^pl /,( r,v)rf 3 


(3) 


is the individual-pressure tensor for the species and f s ( r, v) is 
the velocity distribution function. The superscript in P/ 
indicates that the individual-pressure tensor is defined in a 
frame moving at the bulk flow velocity of species 5. The above 
notation is similar to that used by Rossi and Olbert [1970]. In 
(1) - (3) a and P designate Cartesian components, the particle 
mass is n. = n f = n is the density of the singly charged ions 
and of electrons, p v = m s n is the mass density, and p ? V s a V v p 
is the Reynolds stress or dynamic-pressure tensor. For an 
ion-electron plasma we neglect the small difference between 
the frame moving at the ion bulk velocity V, and the proper 
frame which moves at V, the mass weighted ion plus electron 
bulk velocity defined by pV = . This 

approximation gives P/ f) = P, , where P, is the usual pres- 
sure tensor, defined by (3) with V v replaced by V. An impor- 
tant property of the fluid parameters defined above is that they 
all depend only on the particle information contained in 
f s ( r, v) . All ion fluid parameters were calculated by integra- 
tions over the COT velocity distribution functions. 

The COT ion distribution functions were generated by noting 
that / T (r, v) for any one spatial-velocity or phase space box is 
proportional to the total time all ions spent in the desired velocity 
box while the particles also were located in the spatial box of inter- 
est. For each time step along each trajectory, half the time incre- 
ment was placed in the phase space box in which the particle was 
located at the start of the step and half into the phase space box cor- 


= "v‘?.v £ a + (J S XB )« < 5 > 

The left side of (5) depends only on the particle kinetic energy 
terms defined in (1) - (3). The right side of (5) contains all 
information about the fields and about the fact that a plasma is 
composed of charged particles rather than neutral atoms. The 
sum in (5) is the divergence of the individual -pressure tensor 
V • . The momentum equation describes changes in the 

bulk flow velocity that are produced by force imbalances. The 
next velocity moment of the Vlasov equation gives the energy 
equation, which will be investigated in a separate paper. The 
energy equation is used to study changes in the pressure tensor 
as a result of energy and heat fluxes that flow into and out of a 
given box. As noted previously, the energy equation is needed 
to study topics such as the pressure balance inconsistency 
[Erickson and Wolf , 1980; Hau , 1991; Erickson , 1992; Ash- 
our-Abdalla et al 1 994]. 

Ion and electron equations (5) are added to get a plasma 
momentum equation. When adding ion and electron terms it is 
assumed that the magnetotail is charge neutral and that terms of 
order m e /m i can be neglected relative to unity, e.g. the mass den- 
sity is p = m j n i + m e n e = m- t n . Electrons in the tail are observed 
to have about 1/7 the energy of ions [Baumjohann et al., 1989], so 
that the bulk drift velocities satisfy V g < V. in a quiet time approxi- 
mately 1 -R e thick tail even though the total particle velocities sat- 
isfy v e » Vj- (LK96). Electron drift currents can be more important 
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in thinner tails. The above approximations show that electron 
effects are significant only in the P and j terms of (5) yielding 

3V„ _L dP n „ 

p « + p(V .V ) F a+ X ^ = ( j> < B) a W 

dt p= 1 

where j = j f + j, ; and P pot = P e p „ + P-, pa . As noted previ- 
ously, with the above approximations the ion frame is the 
same as the proper frame which moves at velocity V. Electric 
field forces have cancelled in the neutral plasma. The dV a /dt 
factor in (6) is the acceleration of plasma as seen by a satellite, 
since satellites in the magnetotail are essentially stationary. 
The first two terms in (6) are often combined to give 
p dV a /dt , where dV a /dt is the acceleration of a moving ele- 
ment of the plasma. The second term in (6) is nonzero even in 
a steady state model because the convection speed changes as 
each plasma element flows in from -20 R E to -14 R E (LK96). 
It again is seen that the ion plus electron momentum equation, 
which follows directly from the Vlasov equation, relates parti- 
cle terms on the left side of (6) to particle and field terms on 
the right side of (6). 

It sometimes is useful to think of each individual ion group as a 
separate species with s being the group number. There then would 
be many ion species with the same mass and charge. Equations (5) 
and (6) were used with separate groups of 1000 ions plus electrons 
to check the COT analysis. When using (6) with a single group V is 
the bulk velocity of that group, P is the partial pressure associated 
with that group as defined in the group’s frame, and j is the partial 
current carried by particles in the group. Only B involves the com- 
bination of currents from all particles. 

Kinetic tensors T v can be added since they are all defined in the 
Earth frame. It is not appropriate to weight and add individ- 
ual-pressure tensors when combining the separate groups to form a 
consistent plasma because the P^ for each group was defined in 
a different reference frame. The method used to combine groups 
here and in LK96 was to weight and sum the ion group distribution 
functions. This produced a single ion distribution function for the 
consistent plasma. The final plasma fluid parameters then were cal- 
culated by integrations such as those in (1) - (3) using this com- 
bined ion distribution function. 

Ampere’s law can be used when (6) is applied to the full consis- 
tent model. In this case ( l/ji„) VxB - e () dE/dt can be substi- 
tuted for the combined j from all particle groups, where B and E 
always refer to the preselected fields. The magnetic field cross 
product can be expanded yielding the usual magnetic field pressure 
and tension forces. We conclude that it is possible to fully separate 
particle and field terms only in the final consistent model current 
sheet that is produced by adding contributions from all particle 
groups. The COT method then ideally results in balance between 
the particle terms on the left side of (7) 


3T _L 0Po 

J-jJb 


L \*o B'fi 


(7) 


which were calculated using only f s ( r, v) and the terms on 
the right side of (7), which depend only on the preselected 
fields. 


In this work it sometimes is useful to consider the ions and elec- 
trons separately because ions were treated through orbit tracing and 
electrons through the guiding center approximation. It was seen 
above that the electron and ion nq E forces cancel in a neutral 
plasma. Ion and electron E x B drift currents also cancel. To 
examine ions and electrons separately, it may be noted that these 
statements are equivalent to the observation that the portion of the 
j x B force attributable to E x B drift for either species exactly 
cancels the nq E force for that species, provided E • B =0. One 
simple way to see this involves transforming to a frame in which 
the electric field vanishes at one point of interest, so that both 
E x B and nq E are zero. 

2.3, Electron Terms 

All ion quantities in (6) and (7) were calculated directly. Elec- 
trons were assumed to be isotropic and to obey the guiding center 
equations because the electron adiabaticity parameter k is much 
larger than one. The average ion scalar pressure Pj(x f z) = 
(1/3)7>(P,) was evaluated for each box, and a polynomial fit 
was made to the resulting array. The polynomial fit was used in 
P c (x,z) = ( 1/7) P { (x, z) to calculate the electron pressure. 
The total current carried by guiding center particles 

VP i + •*" tX ( B • V) B 
B 

+ X [ (V, ■ V) V,] (8) 

B B 

was used with P e t[ = P cl to evaluate z) . Equation (8) 
can be obtained by combining guiding center gradient, curva- 
ture, E x B , polarization, and all magnetization drift terms 
[ Parker , 1957]. Although it does not affect calculations of the 
j x B force, the electron parallel current was selected to can- 
cel the ion /„ because a steady state model with negligible y 
dependence cannot support region 1 or 2 current systems. This 
choice of j u also was needed to make V j =0. The last term 
in (8) involves polarization drift, and is negligible for elec- 
trons. Using \ c = \ e /nq e these expressions provide all the 
electron parameters needed in (6) and (7). 

2.4. Tests of Force Balance 

Several tests can be made to see if the model is realistic. The 
electromagnetic force on a unit volume of plasma is written either 
as 



or, when E is uniform, as the right side of (7). The smooth 
solid lines in Figure 1 show the three Cartesian components of 
F EM evaluated using (9) and based solely on the preselected 
fields. The dashed lines are components of j x B calculated 
using the preselected field and currents carried by the COT 
model ions and electrons. The y component of j (jc, z) pro- 
duces the x and z components of j x B . It therefore is not sur- 
prising that the agreement between these components of 
F EM and j x B in Figure 1 are as good as the typical agree- 
ment between the j (a\ z) carried by particles in the model 
and the goal j y (a, z) that is needed to consistently generate 
the preselected magnetic field (LK96). 
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Figure 1 . Each of the 18 panels shows a component of the electromagnetic (9) and j x B forces for a 1 -/? £ wide x 
box as a function of z. A comparison of (6) and (7) shows that these two forces should be equal in a consistent model. 


The v component of j x B in Figure 1 is produced only by j x 
and y,. The least squares fitting routine did not consider j x or j z 
when selecting groups and calculating weighting factors. The parti- 
cle to field energy density ratio or plasma p is near 100 at the equa- 
tor, where ions follow nonguiding center orbits. It is reassuring to 
see that even here the COT technique produces a model with rea- 
sonable balance between all components of the particle and field 
forces. 

Figure 2 separates the electron and ion contributions to the 
j x B force. The x and z components of force from the two species 
are qualitatively similar, but ions clearly dominate. In contrast, the 
middle row of panels in Figure 2 shows that the y components of 
the electron and ion forces both are large but opposite in direction 
so they nearly cancel when added together. The electron contribu- 
tion is based on the smooth polynomial fit to pressures. The jagged 
ion contribution was calculated using the full ion distribution func- 
tion. Noting that different scales are used for the three rows in Fig- 
ure 1, it is seen that the y component of the total j x B force 
matches the desired F^ M = 0 at least as well as the x and z compo- 
nents match the desired F X M and F ^ M . If a model with signifi- 
cant y dependence in the region of interest had been used, the y 
component of the j x B and the electromagnetic field forces would 
be nonzero. 

The reason why ions and electrons in our model have such large 
y components of the j x B force in the sense shown (Figure 2) is 
that the plasma has a strong non-field-aligned earthward compo- 
nent of the bulk velocity even at z = 0 (LK96). An average earth- 
ward drift has frequently been measured by satellites located near 
the neutral sheet [Baumjohann et at. , 1989; Huang and Frank , 
1994]. This cancellation of electron and ion y forces is just what 
one would expect if the earthward bulk velocities are primarily 
E x B drift. 


A real test of the balance between particle and electromagnetic 
forces involves comparing the right side of (7) or (9) with the left 
side of (7). The solid lines in Figure 3 again are F , which 
depends only on the preselected fields. The dashed lines in Figure 3 
are the ion plus electron components of the particle forces, which 
depend only on the distribution functions. The calculated ion pres- 
sure force tends to be a little more jagged than the j x B terms in 
Figures 1 and 2 because the dominant contribution comes from tak- 
ing numerical derivatives of the ion pressure tensor. Even so, errors 
between particle and field terms in Figure 3 are comparable to the 
errors associated with the selection of j (x, z) as shown in Figure 
L This result is important for future work because any study using 
the energy equation depends heavily upon being able to evaluate 
the pressure tensor and other tensors with sufficient accuracy so 
that reliable derivatives can be calculated. 

Electron and ion contributions to each term on the left side of 
(7) have been examined separately. Naturally d/dt is zero in the 
steady state model. The ion contribution to the second term in (7) is 
only about 1% of the following terms for a quiet time model, and 
the electron contribution is even smaller. As discussed later, more 
rapid acceleration is needed to produce injection [DeForest and 
Mcllwain , 1971] and bursty bulk flow [Angelopoulos et al. , 1992] 
events, so the first two terms are significant at times. Since electron 
pressures have been selected to be 1/7 of the ion pressures, it is the 
divergence of the ion pressure tensor that is dominant on the left 
side of (7) and in the dashed curves of Figure 3. 

3. Substorm Effects 

Force concepts are used in this section to study substorms. 
Emphasis is placed on the unbalanced forces that are generated 
when current is disrupted. Substorm studies are complex because 
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Figure 2. Electron and ion contributions to the j x B force are compared in the format that was used in Figure 1 . 


many processes are interconnected. For example, plasma in the 
near-Earth injection region responds to changes in the fields that 
are produced by currents carried throughout the plasma sheet. 
Events such as substorm onsets, injections, and bursty bulk flows 
can be analyzed either in terms of instantaneous local forces and 
the associated plasma acceleration or from the perspective of cur- 
rent disruption and the evolution of polarization electric fields. The 
current system viewpoint is inherently global. It is hard to analyze 
the local substorm effects measured by a specific satellite using the 
current picture without introducing a model of the entire electrical 
circuit. 

Although it is necessary to discuss currents in order to calculate 
forces, some satellite observations can be interpreted more directly 
from the perspective of forces. The j x B and V • P terms in (6) 
give the forces that accelerate a unit volume of moving plasma. 
The local j and B can be measured by particle detectors and mag- 
netometers, respectively, so that changes in the electromagnetic 
force term are directly observable. It is suggested here that the ini- 
tial force imbalance is likely to be produced by a change in this 
electromagnetic term. This suggestion can be tested qualitatively 
by comparing sudden changes in the local plasma flow velocity to 
changes in the measured j x B force. The pressure force cannot be 
measured easily by a single satellite. However, it is likely that the 
principal changes in V - P follow the initial acceleration, taking 
place only after the plasma has moved and therefore become com- 
pressed and heated. One feature of the COT analysis is that fluid 
parameters are evaluated independently in each spatial box. It 
therefore is natural to consider plasma acceleration both in satellite 
data and in the COT analysis from this local perspective. 

3.1. Large-Scale Forces 

Siscoe [1966] calculated the force exerted by the Earth on the 
entire magnetotail by evaluating the x gradient near the Earth of the 
fringing field produced by tail currents. Figure 4 is a sketch sum- 


marizing part of this basic work. The Earth is represented by a 
strong bar magnet with the south pole in the northern hemisphere. 
Each lobe of the magnetotail is treated as a D-shaped solenoid with 
current flowing in the plasma sheet and magnetopause. To illustrate 
forces within the magnetosphere, each lobe was treated as an 
object and replaced by an equivalent bar magnet in Figure 4. Con- 
sidering only the three magnets, it is easy to see that each lobe 
magnet is attracted to the nearest pole of the Earth magnet. This 
illustrates the x component of the electromagnetic force within the 
magnetosphere that keeps the magnetotail attached to the Earth in 
the presence of the tailward force exerted by the solar wind. Most 
of this earthward force is more commonly described in terms of 
field line tension. Attraction of the lobe magnets to each other rep- 
resents the inward electromagnetic force in the z direction that bal- 
ances the outward plasma sheet particle pressure and keeps the 
high p plasma sheet confined near z = 0. An alternate description of 
z forces is that the inward magnetic pressure exerted by the lobes 
balances the outward particle pressure exerted by the plasma sheet. 
The outward magnetic pressure exerted by the lobes on the magne- 
topause is balanced by the inward force exerted by the solar wind. 

Since the Earth dipole field is fixed, any change in solar wind 
stress on the magnetotail requires a change in the current system 
and therefore a change in the strength and/or location of the equiv- 
alent lobe magnets. For example K y , the sheet current density 
integrated through the thickness of the current sheet, increases and 
moves earthward during growth phase, implying that the solar 
wind exerts more force on the magnetotail [Siscoe and Cummings , 
1969]. This is associated with an increase in the total magnetic flux 
in the tail and the resulting increases in tail flaring and x momen- 
tum loss by the solar wind. The COT analysis only modeled a 
small region of the magnetotail so could not be used to calculate 
the total force on the entire magnetotail. Figures 1 and 3 show that 
the z forces are about twice as large as the x forces in our region of 
interest. 
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Figure 3. The two sides of (7) are compared. The electromagnetic force depends only upon the preselected fields and 
the particle force depends only upon the /( r, v) calculated during the COT analysis. 


3.2. Substorm Forces and Plasma Acceleration 

Most steady state studies of the quiet magnetotail neglect the net 
unbalanced force needed to accelerate plasma even though the 
plasma p is near 100 at the center of the current sheet. A compari- 
son of terms in (6) shows the reason this omission is justified. A 
plasma element is accelerated earthward by the difference between 
the earthward force on a unit volume given by the x component of 
j x B or j B z in the 2-D model and the tailward force provided by 
the x component of V • P . Plots of the first two terms in (6) were 
not shown separately because the first is zero and the second is 
much smaller than the j x B and V • P terms in the steady state 
COT model. LK96 showed that V x near z = 0 decreased by 50 km/s 
as one moved earthward by 5 R E , yielding an average acceleration 
of 0.2 km/s 2 . The net force needed for this weak quiet time acceler- 
ation in the COT model is 1 0 -19 N/m 3 , which represents only a 1 % 
difference between the j x B and V * P terms shown in Figure 3. 

If force imbalance and the associated acceleration always was 
as weak as in the quiet time COT model, it would take more than 
half an hour for an initially slowly moving plasma element to reach 
a speed of 400 km/s. During this time the plasma would move 
more than 60 R £ . It is not known how long it takes for plasma in 
rapid flow regions to accelerate. Direct measurements of the 
instantaneous acceleration of a moving plasma element would be 
difficult to make even with an array of satellites. The mere obser- 
vation of high speed flow also does not require a large local accel- 
eration since no net force is needed to maintain even rapid steady 
flow. As an example, the net force on a rapidly moving mature 
plasma bubble can be zero [Chen and Wolf * 1993; Sergeev et a/., 
1996c]. However, injection events require strong plasma accelera- 
tion throughout the midnight region and any observed rapidly flow- 
ing plasma must have been accelerated somewhere. Accelerations 


and force imbalances well above those in the quiet time COT 
model therefore must sometimes be present in the magnetotail. 
Borovsky et al. [1997] emphasized that the observed instantaneous 
plasma velocities often fluctuate very rapidly in both magnitude 
and direction even during quiet times. This observation implies 
that large force imbalances may be common and therefore impor- 
tant to any study of the magnetotail. 

As an example of the very rapid acceleration that may some- 
times be present, it would take 40 s for an initially stationary 
plasma element to reach 400 km/s if the average acceleration were 
10 km/s 2 . The plasma would move 1.3 R E during this time. The 
force required for such an acceleration of the COT plasma is F x = 
5 x 10“ 18 N/m 3 , which is a large fraction of the j x B and V • P 
terms in Figure 3. We conclude that significant differences between 
V • P and j x B are needed if such large acceleration events take 
place either during substorms or during normal conditions in the 
middle magnetotail. 
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3.3. Substorm Currents 

In a model with little y dependence, changes in K y extend from 
flank to flank and close over the northern and southern magneto- 
pause. A sketch of the northern lobe current system is shown in 
Figure 5a. This is the usual picture of growth phase. Observations 
suggest that K y initially is reduced at substorm onset only over a 
limited (x, y) region, with the current loop closed by field aligned 
current to the ionosphere (Figure 5c). The forces associated with 
this mode of current disruption are considered here. 

The previous discussion shows that the earthward component of 
either j x B or V • P must change substantially to produce the 
plasma acceleration seen during substorms. Substorm onset usually 
is associated with a large change in magnetotail currents. This cur- 
rent change will produce changes in B throughout much of the 
plasma sheet. As suggested above, it therefore appears that changes 
in the j x B force should be examined in order to understand the 
cause of the initial substorm plasma acceleration. 

Magnetic field changes have been used to estimate the location 
of a satellite with respect to the current disruption region [Lopez et 
al ., 1988; Ohtani et al., 1992; Lui, 1996]. A reduction of K y in a 
given (x, y ) region results in a decrease of B x in the outer current 
sheet. A more important effect of a localized reduction of K is the 
resulting increase of B z everywhere inside the diversion loop as 
shown in Figure 5d, The increase of 5, in the inner magnetotail is a 
reliable observational signature of substorm onset. The assumption 
that K decreases only in a localized region implies that K y ini- 
tially remains constant while B 2 increases earthward of the disrup- 
tion region. Although the shape of magnetic field lines will change 
everywhere, resulting in a change in the particle drift velocities, the 
average over an extended region cannot change greatly in a 
tail-like structure. A certain average tail lobe B x is needed so the 
outward lobe magnetic pressure balances the inward solar wind 
pressure at the magnetopause. A certain average AT is in turn 
needed to generate this lobe B x . A large increase in B z combined 
with a relatively steady K earthward of the disruption region 
therefore produces an increase of K y B z , the earthward force per 
unit area of the current sheet. At the same time B decreases, the 

z 

earthward force decreases, and field tines become more stretched 
outside a diversion loop (Figure 5d). 

a) b) 



f i l 


c) 




Figure 5. (a) Growth phase currents in the northern tail lobe, (b) 
view looking down on the equatorial plane of the growth phase 
currents and the resulting equatorial magnetic field, (c) perturba- 
tion currents that are added to those in Figure 5a at substorm onset, 
(d) view looking down on the equatorial plane of the substorm cur- 
rent loop and the resulting perturbation fields. 


Pulkkinen et al. [1991] showed quantitatively how much the 
equatorial end points of field lines with fixed ionospheric foot 
points both inside and outside a current loop move, and Kaufmann 
[1987] showed how much foot points of field lines that connect to 
fixed tail locations move when current is diverted in a model mag- 
netosphere. These two studies added extra current systems to exist- 
ing magnetotail models. The nearly 2-D COT model that was used 
here is adequate to estimate the magnitude of the force changes 
needed to produce the observed substorm plasma acceleration near 
the center of a diversion loop. However, 3-D simulations or other 
time dependent techniques are needed for a comprehensive analy- 
sis of the current diversion process and of substorm injection 
events [Sdnchez etal 1993]. 

A reduction of K through the formation of a diversion loop 
generates the important increase in B , more efficiently than would 
a decrease of K from flank to flank with closure on the magneto- 
pause. Figure 5d shows that both the reduction of K and the asso- 
ciated field aligned currents enhance B z everywhere inside a 
current diversion loop. In contrast, the east-west currents in the 
current sheet and magnetopause are in opposite directions when K 
extends from flank to flank (Figure 5b). If this typical growth phase 
current system simply decreases, then A K y in the current sheet 
increases B „ and the change in the magnetopause part of the current 
loop decreases B z at the equator earthward of the current loop. The 
current sheet effect dominates at z = 0 so that a small net increase 
of and the associated earthward force are produced in the earth- 
ward region at the equator by the fringing field of the full loop as 
the current in Figures 5a and 5b decreases. 

A very simple current carrying wire model can make the above 
conclusions more quantitative. Earlier studies [ Kaufmann , 1987; 
Kaufmann and Larson, 1989 \ Kaufmann etal. 1990; 1993] showed 
that reasonable estimates of magnetic field perturbations could be 
made with simple wire models provided the observer remained 
about one scale length away from the wires. The example used 
here models current disruption well out in the tail, but illustrates in 
general how sensitively A B. depends on the current closure path. 
First, the Biot-Savart law was used with one of the loops shown in 
Figure 5a to calculate A B when the current is reduced from flank 
to flank and the loop closes around the semicircular magnetopause. 
The resulting A B z is (32” l/2 ) (\i 0 AI/nR ) s A B 0 as seen by an 
observer located at midnight in the equatorial plane at a distance R 
either earthward or tail ward of the localized disruption loop. In this 
expression A / is the current added to the wire loop and R is the 
radius of the tail. Current disruption increases B z by A B 0 at a dis- 
tance R earthward of the loop and decreases B by A B 0 at a dis- 
tance R tail ward of the loop. 

A simplified version of Figure 5d was used to make the compar- 
ison with current diversion to the ionosphere. A flat rectangular 
loop in the equatorial plane was used for the simple estimate 
instead of the circuit in Figure 5d which has current flowing along 
curved field lines. The rectangular model appears to be reasonable 
if disruption takes place in a thin plasma sheet. The contribution 
from the ionospheric leg of the rectangle was neglected because it 
was assumed to be far from the observer. The resulting increase in 
A B. is ( 1/2 + 2 -I/2 ) (\i o AI/kR ) = 6.8 A B a when current dis- 
ruption extends from flank to flank but the loop closes in the equa- 
torial plane and the midnight observer is a distance R earthward of 
the cross-tail disruption wire. If the current disruption rectangle is 
only half this wide, extending a distance R/2 from midnight toward 
each flank, then the midnight observer located a distance R earth- 
ward of the cross-tail leg sees A B z increase by (1 + 1.25 1/2 ) 
(\i () AI/nR ) = 12 AB n . For an observer at midnight a distance R 
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tailward of the cross-tail leg, decreases by 1.2 A B () and 0.7 
A B a using the rectangles of half width R and Rf 2, respectively. We 
conclude that current diversion to the ionosphere can be an order of 
magnitude more efficient in producing A B z and the associated 
earthward force on plasma inside the loop than a simple reduction 
of the growth phase current system. There is a much weaker depen- 
dence on the mode of disruption for plasma located tailward of the 
disruption region. 

In summary, the principal magnetic signature of local dipolar- 
ization, the observation of a sudden increase in implies that a 
satellite is inside the principal diversion loop [Lopez et al . , 1988; 
Ohtani et al., 1992; Lui, 1996]. This field change results in a net 
earthward force on a unit volume if j y (x, z) remains nearly con- 
stant at the observation point, i.e., if current is disrupted in a local- 
ized region tailward of the observer. On the other hand, a satellite 
that is just outside the diversion loop will see the magnetic field 
rapidly become more stretched. This is equivalent to a sudden 
enhancement of the growth phase perturbation, as noted by Ohtani 
et al. [1992]. Although the resulting acceleration of a moving ele- 
ment of plasma is difficult to measure, it should be possible to mea- 
sure the change in j x B and therefore to calculate the expected 
initial acceleration at substorm onset. 

3.4. Pseudoexpansions and Multiple Expansions 

Several proposed disruption mechanisms require a thin current 
sheet. Such mechanisms can create new diversion events on the 
newly stretched field lines tailward of and at the sides of an original 
diversion loop. This sequence is similar to that seen during multi- 
ple expansion substorms [Rostoker, 1996; Sergeev et al ., 1996b]. 
Figure 6 is a sketch of a possible multiple expansion process. It is 
assumed that the current increment labeled Aj, is the first to be 



Figure 6. Growth of a multiple expansion substorm. Initiation of 
the first substorm current loop Aj, causes a thinning of the plasma 
sheet everywhere outside this first loop. If nearby regions become 
thin enough, additional substorm loops Aj 2 and Aj 3 may be cre- 
ated. 


produced when an appropriate localized region of the tail becomes 
sufficiently thin. This current change causes regions adjacent to the 
initial diversion loop to become thinner (Figure 5d). A disruption 
process that requires a thin current sheet in an appropriate tail 
region then can generate the Aj 2 and Aj 3 diversion loops in Fig- 
ure 6, thereby extending the disruption region. A full substorm may 
require the formation of a sequence of diversion loops. A pseudo- 
expansion would result if only a small region of the tail becomes 
thin enough in a region that can cause disruption. 

It was previously noted [ Kaufmann , 1987] that nearly all ions 
remain in the magnetotail or outer magnetosphere during the brief 
period of substorm onset. A substantial fraction of the magnetotail 
electrons can be exchanged through field aligned currents with 
electrons in the ionosphere. However, if the total number of tail 
ions does not change suddenly at onset, then quasineutrality 
requires that the total number of electrons does not change either. 
Ions carry most of the cross-tail current in the COT current sheet 
model. Models containing very thin current sheets just before sub- 
storm onset [Lee et al., 1995; Pritchett and Coroniti , 1995; Ma et 
al., 1995; Birn etaL, 1996] can have much of the cross-tail current 
carried by electrons. The above comments suggest that a study of 
substorm onset should consider mechanisms that reduce the 
cross-tail current carried by either species tailward of the injection 
region without decreasing the total number of particles in the mag- 
netotail. Kaufmann et al. [this issue] describe such processes 
involving nonguiding center motion in a thin current sheet. Causes 
of the initial plasma sheet thinning also are discussed by Kaufmann 
et al. [this issue]. 

4. Discussion and Summary 

A consistent orbit tracing (COT) method was used to find the 
particles needed to produce a quiet time model magnetotail. The 
model has little y dependence and a characteristic thickness scale 
length of 0.7 R E in the -20 R E < x < -14 /? £ , 0 < Izl < 2 /? £ region 
of interest. Orbit tracing started with an isotropic Maxwellian pop- 
ulation of protons. Electrons were added using the guiding center 
approximations. Several similar nearly consistent current sheets 
were created in our earlier study (LK96) starting with monoener- 
getic protons at energies between 1 .5 and 1 5 ke V. 

4.1. Uniqueness 

Appendix A addresses some questions concerning the unique- 
ness of COT model results. It is concluded that the 2-D tail models 
used here are not uniquely generated by a specific distribution 
function. However, it was found that groups of nonguiding center 
ions following each type of orbit are required. The selection of par- 
ticle energies for various runs was based on satellite observations. 
If only 5-keV protons that mirror at low altitude (Speiser orbits) are 
present, then a current sheet is formed with a thickness of approxi- 
mately z = 0.25 R e . The energy dependent parameter z Q = 

mv/[qB x (z 0 )] is the point at which the particle’s distance from 

the equatorial plane is equal to the z component of the particle s 
local gyroradius. In the above expression m, v, and q are the parti- 
cle mass, velocity, and charge, and B* y = B 2 + B 2 . The analysis 
showed that the addition of trapped particles can broaden the cur- 
rent sheet to a thickness of about 2z„ . Trapped particles meander 
back and forth across z = 0 but never enter the region in which 
orbits spiral around field lines and particles are magnetically mir- 
rored. The 2z„ scale still is thinner than is required for the model 
in the -20 R E < x < -14 R £ region of interest. Particles on cucum- 
ber orbits that mirror throughout the current sheet also had to be 
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included to produce the consistent current sheet. We were not able 
to generate a consistent quiet time model current sheet unless parti- 
cles following all three orbit types were included. Therefore, 
although a unique /(r, v) is not obtained, the mixture of particle 
reflection or mirror points is constrained if the particles are to carry 
the j y ( x , z) that is needed to produce the model magnetic field. 
The more difficult question of uniqueness in a 3-D model with 
strong y dependence is not addressed here. 

Appendix B addresses another aspect of the question of unique- 
ness. The problem considered is how charge neutrality can be 
assured all along a field line on which electrons follow guiding 
center orbits and ions follow nonguiding center orbits. The COT 
method used isotropic electrons at the equator plus a parallel elec- 
tric field to produce a neutral plasma. This method of introducing 
an E |( to maintain charge neutrality had been used in several early 
studies [Alfvtn and Fdlthammar , 1963, p. 162; Lemaire and 
Scherer , 1971; Chiu and Schulz , 1978]. Appendix B presents an 
alternate technique. An explicit method is developed to calculate 
the equatorial electron pitch angle distribution that will assure 
charge neutrality all along a field line with E % = 0 for a given varia- 
tion of the ion density along the field line. We conclude that it is 
possible to pick a variety of electron distributions and electric 
fields in a charge neutral current sheet, so that neither E n nor the 
electron distribution function is uniquely specified by the selection 
of a magnetic field model. 

4.2. Force Balance 

The set of 120 distribution functions generated by the COT 
analysis was treated much as if one had data from a uniformly 
spaced array of 120 satellites. Evaluating the terms involving spa- 
tial derivatives in the momentum equation (6) requires a knowl- 
edge of fluid parameters throughout the region of interest. It 
therefore was easy to examine all the force and acceleration terms 
in the model calculations. In contrast, evaluating V • P and 
( V • V) V are difficult when dealing with satellite data. The j x B 
and d/dt terms in (6) are easier to measure because they require 
data at only one location. 

Figures 1 to 3 were used to examine force balance. Since j car- 
ried by the particles was considered when combining plasma 
groups, it was not surprising that the resulting x and z components 
of the j x B force were almost equal to the electromagnetic force 
required to create a consistent magnetotail (Figure 1). Even though 
j x and j z were not used in the fitting procedure, it was found that the 
resulting y component of the j x B force also agreed with the 
requirements for consistency. Figure 2 showed that ion forces dom- 
inated the x and z components of j x B , but electron and ion terms 
almost cancelled each other in the y direction in the nearly 2-D 
modeling region. This cancellation was understood as primarily an 
effect of earthward drift which is associated with E . Figure 3 
showed that the V • P force balanced the electromagnetic force to 
a good approximation in the model. This result demonstrated that 
the full pressure tensor is accurately calculated during a COT anal- 
ysis. Inertial forces associated with the bulk motion of plasma par- 
ticles were much smaller than the pressure forces in the quiet time 
model. 

It is useful to compare the COT results to previous kinetic stud- 
ies. The companion paper examines individual elements of the 
pressure tensor. Here we only consider the net forces. Burkhart et 
al. [1992] used a 1-D simulation. The ion source in this study was 
located far from z = 0, so there were few trapped particles. Forces 
were found to be well balanced in thin current sheets with a small 
adiabaticity parameter k. Force balanced solutions could not be 


obtained in thin current sheets when the average ion k approached 
one. Problems were encountered in balancing the z component of 
particle and field forces. Our analysis required large contributions 
from ions on trapped and cucumber orbits in the thick quiet time 
current sheets associated with large k. Burkhart et al. [1992] noted 
that the force balance problem in thin large k current sheets may be 
associated with the nature of the 1-D model with few trapped ions. 
This issue is examined in more detail in the companion paper. 

Ashour-Abdalla et al. [1993, 1994] used large-scale kinetic 
(LSK) orbit tracing to model all the magnetotail earthward of a dis- 
tant neutral line. The magnetic field model was 2-D and based on 
T89. A drifting Maxwellian mantle source was used and orbits 
were traced for such a long time that Speiser, cucumber, and 
trapped orbits all became well populated before particles reached 
the middle magnetotail. Current sheets in this model tended to be 
thinner than is needed to produce the original T89 field. The thick- 
ness of the principal model current sheet was comparable to the 
thickness of the region in which trapped ions meander. 

Pressure balance was examined in the LSK model. Forces were 
found to be reasonably well balanced, especially in almost 1-D 
regions of the tail. The force balance test used was a combination 
of the marginal firehose stability condition beyond the edge of the 
current sheet, which involves x forces, and the need for z forces to 
balance the lobe pressure. The guiding center approximations were 
well satisfied in the outer current sheet, so only the P\ \/P ± ratio 
was needed in this test. Kaufmann et al ., [this issue] show that 
other elements of P are needed for a detailed study of forces very 
close to z = 0. 

4.3. Substorm Effects 

Substorm onset was studied as an application of the force bal- 
ance analysis. The COT results pointed out that the average plasma 
acceleration in a quiet current sheet requires only small differences 
between the earthward j x B force and the tailward V • P force 
but that bursty bulk flow and substorm injection events require 
more substantial force imbalances. Substorm onset involves a 
localized reduction of the integrated cross tail sheet current density 
K y near or somewhat tailward of that portion of the plasma sheet 
that begins to collapse. The reduction of K must take place with- 
out significantly reducing the total number of particles in the tail 
during the brief onset period. Nonguiding center mechanisms that 
can do this are discussed by Kaufmann et al. [this issue]. The 
increase of B z and the associated increase in the earthward K y B z 
force were treated as immediate causes of the initial plasma accel- 
eration. It was found that diversion of cross-tail current along field 
lines to the ionosphere, which involves reduced K in a region that 
is limited in both the x and y directions, is much more efficient at 
producing an increase in B z and the net earthward force than is a 
reduction of K y that extends from flank to flank and then closes on 
the magnetopause. 

Cross-tail current disruption and the associated current diver- 
sion loop cause stretching of field lines to the east, west, and tail- 
ward of the diversion loop. If disruption is produced only in very 
thin current sheets Kaufmann et al. [this issue], then new diversion 
loops can be generated adjacent to the original loop. This mecha- 
nism provides one way to produce multiple expansions or the 
spreading of a substorm. 

Appendix A: Uniqueness 

One feature that the 2-D magnetostatic COT solutions have in 
common with electrostatic BGK [Bernstein et al ., 1957] solutions 
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is that neither method finds a unique distribution function. This 
appendix describes some of the simplest nonunique aspects of the 
analytic magnetic field model that was used. Nonunique aspects of 
the COT ion distribution functions and of the techniques used to 
treat electrons also are described. 

AL1. Harris-like Field Models 

The original 1-D Harris [1962] model and the many 2-D Har- 
ris-like variations [Kan, 1973] start with an assumed /( r, v) at z = 
0. All these models are based on the fact that any function that can 
be written in terms of only constants of the motion is a solution to 
the Vlasov equation. A drifting Maxwellian at z = 0 was used for 
the I-D and 2-D cases because it is easily written as a function of 
the hamiltonian H and the y component of the canonical momen- 
tum P , both of which are constants of the motion. Kan [1973] 
found that an infinite number of solutions for the 2-D magnetic 
field could be found with the assumed /(r, v) , all of which satisfy 
Maxwell’s equations. Three dimensional cases are more difficult 
because P (r, v) is not a constant of the motion. The uniqueness 
of 3-D analytic models is not considered here. 

The Harris-like set of models was derived using some assump- 
tions that are not realistic in the magnetotail. One assumption is 
that a reference frame exists in which E (r) =0 everywhere. The 
assumption that most clearly disagrees with observations is that the 
electron and ion cross-tail bulk velocities V and V jy both are 
independent of r throughout the entire plasma sheet. The observa- 
tions and COT results described in LK96 involved a V iy (r) that 
increases as one moves toward z = 0 and toward the Earth. Finally, 
the Harris-like magnetotail model we used was based on an expan- 
sion in which certain small terms were assumed to be negligible 
[Bimetal., 1975; Zwingmann, 1983]. 

Perhaps the simplest way to see that the 2-D Harris-like fields 
can be generated by different distribution functions involves the 
dimensionless variables that are used to solve Ampere’s law [Kan, 
1973]. The distance r is normalized using the factor h defined by 
2\i 0 h 2 = T/n a q 2 V 2 where V = V jy -V ey ,T= T f + T e is a ther- 
mal energy, and n n is a density normalization constant. The vector 
potential is normalized using a = 2h (Snn a T) 1/2 . Two vector 
potentials Aj (jc, z) and A 2 ( x , z) therefore are equal at the same 
unnormalized (x, z) point if n o\ T \ ~ n o2^2 and n ol^[ ~ 
n 0 2 V 2 so that h } - h 2 and a f - a T For example, one could gener- 
ate the same magnetic field by doubling the thermal energy and 
drift velocity of all ions and electrons while reducing the density by 
half since this does not alter either r or A ( r) . More important, the 
same A (jc, z) can be produced by adding n ol / 2 of the particles 
discussed above with T lf V x to n 0 2 /2 of a group with T 2 , V 2 * 
Similarly, this same A ( x , z) can be produced by adding different 
amounts of many groups of drifting Maxwellians, each having the 
same T/V ratio. A wide variety of distribution functions can be 
generated by combining various amounts of several drifting Max- 
wellians. Each /( r, v) produced in this manner generates the 
same 2-D Harris-like B (r) . 

A1.2. COT Results 

Rather than starting with a preselected /( r, v) and solving for 
B (r) , the COT procedure started with a preselected B (r) and 
searched for /( r, v) . A uniform cross-tail electric field was 
included for the case studied here. LK96 showed that satisfactory 
approximate solutions could be found for a given B (r) by start- 
ing with monoenergetic ions, even though the original Harris-like 
model was derived using a drifting Maxwellian at z = 0. LK96 also 
found that particles with different energies could produce models 


that are nearly equal in approaching the goal of consistency. Fig- 
ures 1 to 3 show that results of comparable consistency are 
obtained when the COT method is started using Maxwellian ions. 

Although the COT method produced reasonable plasma sheets 
for nearly all 1-D and 2-D magnetic field models we have used to 
date, there are some models that cannot be produced using particles 
from a preselected parent population. For example, thin current 
sheets cannot be generated using ions that are so energetic that they 
cannot be deflected significantly by the magnetic field within the 
current sheet. As an extreme example, a model with a discontinu- 
ous B could not be generated by any finite energy particles because 
the model would require an infinitely thin current sheet. The thin- 
nest current sheets that can be produced by particles with a given 
energy have a thickness of approximately z Q . Kaufmann et al. [this 
issue] also show that we could not create thin current sheets using 
only particles with energies that produce the most chaotic orbits. 
Only limited tests have been carried out to date in an attempt to 
find particles that generate 3-D model fields with strong y depen- 
dence. An extensive study of 3-D cases will be more challenging 
than the quiet time 2-D example used here. However, current in the 
actual magnetotail is carried by real ions and electrons. As a result, 
it should be possible to find a distribution of particles that will gen- 
erate a sufficiently realistic model of the magnetotail. 

With regard to the possibility of finding radically different parti- 
cle distributions that produce the same magnetotail field, it was 
noted above that similar fields could be produced using a variety of 
quite different energy distributions. However, once the energy dis- 
tribution is selected to agree with observations, then the COT stud- 
ies found that the distribution of particle orbit types and reflection 
points was constrained. For example, we have traced many groups 
dominated by trapped particles and many groups dominated by 
Speiser particles. In every case the trapped particles show the char- 
acteristic negative j near z = 0 and positive j y in the z () < Izl < 2 z () 
region. Similarly, Speiser groups uniformly carry positive j y 
peaked at Izl < z fl . There always will be a range of particle distribu- 
tions that produce almost equally good current sheets using the 
COT method. For example, it is well known [e.g., Stix, 1962, p. 
114] that a Vlasov plasma can be approximated as closely as 
desired by combining many beams that follow free-streaming or 
unperturbed particle orbits. Highly structured distribution func- 
tions containing multiple beams have been proposed and observed 
in some parts of the plasma sheet [Frank et al, 1994; Ash- 
our-Abdalla et al ., 1996]. It is unlikely that the COT method will 
be able to select one specific highly structured distribution function 
from a group of structured distribution functions that produce simi- 
lar fluid parameters. 

A 1.3. Electron Assumptions 

Only ion orbits were traced in the COT analysis because ther- 
mal electrons obey the guiding center approximations in the model 
used here. If the magnetic moments of both ions and electrons were 
conserved, then steady state charge neutrality could be assured 
everywhere along a field line by selecting £ H = 0 and the same pitch 
angle distribution for ions and electrons at the equator. However, 
the mirroring process will not assure plasma neutrality along a 
magnetotail field line when electrons follow guiding center orbits 
and ions follow nonguiding center orbits. Plasma was kept neutral 
along field lines in the COT calculations by using the Boltzmann 
relation to determine the required parallel electric field. Electrons 
usually were assumed to be isotropic and Maxwellian at z = 0 
[Kaufmann and Lu , 1993]. The resulting small parallel electric 
field primarily modified the electron density to produce neutrality. 
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The density of the more energetic ions was almost unaffected by 
this E r A small relatively unimportant perpendicular electric field 
in the x-z plane also was needed to keep VxE = 0 in this steady 
state 2-D tail model. 

The above electron assumptions are not the only way to produce 
a neutral plasma. Appendix B describes a method of finding a 
nonisotropic pitch angle distribution for guiding center electrons 
that will assure plasma neutrality through simple mirroring effects 
with no parallel electric field. This technique was not applied in the 
COT runs carried out to date because the ion density was kept only 
in the rectangular region of interest near the equator. The ion den- 
sity is needed all the way down to the Earth to use the technique in 
Appendix B. This appendix does show that the parallel electric 
field used in the analysis is not unique. A wide variety of electron 
pitch angle distributions could be combined with various E n 
choices to maintain charge neutrality. However, very strong pitch 
angle dependence would be needed to cause electron densities to 
change as rapidly as do the ion densities near z = 1 R E in Figure 
10a of LK96 if this density change is to be produced by the mirror 
effect alone. 


The purpose of this appendix is to invert (B5) to get the 
unknown pitch angle distribution (0 t ) from n (5) , the known 
density distribution along a field line. The solution can be 
described most clearly by changing variables several times. First 
setting £ = sin 2 0, gives 


+ (£) = f F|ISin < B6 ) 

Jo ($-0 

Then multiplying both sides by (o-^) W2 dt, where o is a 
new variable and integrating from 0 to o gives 
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Appendix B: Charge Neutrality 

This appendix shows how to find the pitch angle distribution of 
guiding center electrons that is needed for a plasma to be charge 
neutral everywhere along a magnetic field line. It is assumed that 
£■„ = 0 and that the density of nonguiding center ions is known all 
the way from a reference point at which B = B { down to the Earth, 
where n. goes to zero. The reference point can be the equator in the 
middle magnetotail because electrons follow guiding center orbits 
all along the field lines. The variable 5 = B^/B is used to specify 
the location along a field line. The subscript 1 refers to the refer- 
ence point where 5=1- The electron magnetic moment must be 
conserved so that 


where the order of integration has been reversed in the final 
form. This has been done because the d^ integral is just equal 
to 7 i, as can be seen by changing variables from 5 to t = 
(5“C)/(° _ C) wh ere C is fixed for this integral. Equation 
(B7) then becomes 

f = ”f fi[sin _l (C l/2 )]^ (B8) 

Jo < o - 5 ) Jo 

The inversion is now completed by taking d/do of both sides 
of (B8). A convenient form for the left hand side is obtained 
by changing variables from 5 to P = o- 5 before taking the 
derivative. This gives 


sin 2 0/£ as sin 2 0 ] /5 1 (Bl) 

where 0 is the pitch angle. The region of interest is 5 ^ 1 • The 
pitch angle distribution at the point 5 is defined as 

F^ (0) = f°2JCv 2 /(S, v, 0) dv (B2) 
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where the velocity distribution function /(5, v, 0) is assumed 
to be independent of the phase angle around B. The particle 
number density is 


where (0) =0, the primes indicate differentiation, and the 
last form has been transformed back from p to 5- Setting o = 
sin 2 0j gives the desired F, (0|) 


-jc/2 

n (£) = 2 1 Ft (0) sin0d0 (B3) 

Jo 

assuming symmetry about 0 = tc/ 2 . Liouvi lie's theorem 


-5 


» sin 2 8 

Jo 


1 j»(5)+25«'(5) * 


(BIO) 


F k (6) = F, (0|) = F, [siiT l (£ 1/2 sin0)] (B4) 

relates the distribution function at 5 to the distribution func- 
tion at the reference point. Substituting (B4) into (B3) and 
changing variables using (Bl) gives 


A final form without the apparent problems at the limits can 
be obtained by changing variables from t, to a = a -5 * n 
(B9) before setting a = sin 2 0j 
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where the derivative is with respect to the argument 
sin 2 0, - a 2 . Equation (Bll) can be evaluated for any 0, pro- 
vided the particle density n (#,//?) is known for all B> B f . 
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